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1. Report Summary 

Reported here are the results of our numerical investigation into the mechanisms which 
affect the transport and growth processes in Physical Vapor Transport (PVT) crystal 
growth ampoules. The first part of the report consists of a brief summary of the major 
accomplishments and conclusions of our work. The second part consists of two manuscripts, 
submitted to the Journal of Crystal Growth , which provided a detailed description of the 
findings in our investigation. 

2. Project Objectives 

The objectives of the project were to conduct a detailed numerical examination of the 
vapor transport mechanisms which govern the growth of crystals in PVT ampoules. A 
specific focus of the project was to examine the effects of thermal creep and thermal stress 
convection on the crystal growth process. These mechanisms result from temperature 
gradients along the solid surfaces of the system (thermal creep) and within the bulk gas 
(thermal stress). Typically, the convective velocities resulting from thermal creep and 
stress are significantly smaller than those resulting from buoyant convection under similar 
nonisothermal conditions in normal gravity. However, in microgravity conditions the creep 
and stress mechanisms could potentially have a significant effect on vapor transport relative 
to diffusion-limited conditions. 

A second focus of the project, which was developed subsequent to the initial proposal, 
involved the examination of the effects of conduction heat transfer within the solid crystal, 
and deposition of the vapor onto the side wall the ampoule, on crystal growth rates and 
growth uniformity. 


3. Effect of Thermal Creep and Stress on PVT 


Thermal creep is a slip flow of gas over a solid surface which results from a temperature 
gradient tangential to the surface. On a phenomenological basis, the thermal creep velocity 
is given by 


C.,v dT 



dz 


5 

wall 


(i) 


where v is the kinematic viscosity of the gas, z is the coordinate tangential to the wall, 
and C 3 is a dimensionless coefficient , of order unity, which depends on the momentum and 
thermal accommodation of the molecule-wall collisions. In ampoules with nonisothermal 
side walls, thermal creep would result in a recirculating flow within the ampoule. The gas 
would slip along the side walls towards the higher-temperature source, and return along 
the centerline to the lower-temperature crystal. 

To examine the effects of creep-induced recirculation, we first developed a numerical 
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model to predict the vapor growth rates in PVT ampoules. The model consisted of a 
steady-state, two-dimensional axisymmetric representation of a binary component (nutri- 
ent + inert carrier) PVT system. A fully compressible formulation was employed in the 
formulation of the governing continuity, momentum, energy, and species equations. Along 
with thermal creep and natural convection, the model included transport of nutrient by 
Soret diffusion. To facilitate the identification of regimes in which thermal creep had a sig- 
nificant effect on growth rates, a non-dimensional formulation of the governing equations 
was adopted. 

Details of the formulation and calculation results appear in Mackowski et al. (1996b) 
(given in Appendix A). We found that thermal creep will have a significant effect on 
the distribution of mass flux at the crystal interface only for relatively large Schmidt 
number (Sc — i '/D ab, where v and Dab are the kinematic viscosity and binary diffusion 
coefficient) conditions, and relatively large wall temperature gradients adjacent to the 
crystal interface. Such conditions, however, have been realized in a number of ground and 
space-based PVT experiments, in particular those involving systems containing a trace 
concentration of a large molecular mass nutrient in a low molecular mass carrier (such as 
the 3M PVTOS experiments). 

The overall effect of thermal creep convection on crystal growth is qualitiatively sim- 
ilar to that resulting from buoyant convection. Thermal creep creates a recirculating 
flow within the ampoule, with the fluid driven along the ampoule walls from the lower- 
temperature crystal to the higher-temperature source, and returning to the crystal along 
the ampoule centerline. The recirculation augments the deposition rate at the crystal sur- 
faces nearest the centerline, and retards the growth rates near the side wall. Variations in 
mass flux of up to 50% in the presence of thermal creep were identified for realistic PVT 
growth parameters, as opposed to essentially uniform mass flux for creep and buoyancy- 
free conditions for the same operating parameters. 

Even in systems where the absolute difference in the source and crystal temperatures 
is not large, relatively large wall temperature gradients adjacent to the crystal can occur as 
a result of wall heating strategies that are used to control parasitic nucleation on the walls. 
In conjunction with thermal creep, these gradients can lead to gas recirculation patterns 
directly above the crystal interface, which in turn can result in a distribution of flux at 
the interface that is significantly more nonuniform than that predicted in the absence of 
creep. 

We also performed a limited number of calculations that included the convection mech- 
anism of thermal stress. Thermal stress is predicted from the Burnett contributions to the 
fluid stress tensor — which constitute the second-order approximation to the Boltzmann 
equation for small Knudsen flows (the Navier-Stokes equations being the first order ap- 
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proximation). These constitutive relations predict that a non-linear temperature gradient 
in the fluid gas can directly impart a fluid stress, and lead to a bulk convection of the fluid. 

Our predictions indicate that 1) thermal stress convection is only significant for ex- 
tremely high temperature gradient conditions (i.e., ~ 100 K/cm), and 2) the thermal creep 
flows generated from these temperature gradients invariably overwhelm the thermal stress 
flows. Consequently, in regard to crystal growth, we do not believe that conditions exist 
in which thermal stress, by itself, would have an appreciable effect on the growth rates 
and morphology of the crystal. However, the thermal stress flows that we have predicted 
are intriguing in themselves. Under earthbound conditions thermal stress convection in 
nonisothermal gases will be insignificant compared to buoyant convection. The /ig environ- 
ment, on the other hand, offers conditions in which thermal stress convection could become 
a significant if not dominant mode of convection in nonisothermal gases. Our future work 
will focus on thermal stress convection in single-component (i.e., no mass transfer), near- 
continuum gases, with the goal of identifying designs for microgravity-based experiments 
to measure thermal stress convection in nonisothermal gases. 

4. Effect of solid— phase heat transfer and wall deposition 

Previous PVT numerical modelling efforts (including those discussed above) have 
taken the crystal interface to be planar and isothermal. In reality, conduction heat transfer 
through the solid, coupled with the latent heat release of the depositing crystal, can lead to 
nonuniformities in the crystal surface temperature. This, in turn, will affect the saturation 
partial pressure of the nutrient above the crystal and thus alter the arriving flux. In 
addition, parasitic nucleation on the side walls of the ampoule can lead to a growth interface 
that, overall, is non-planar. 

To examine these issues, we modified our numerical model of the PVT system to 
include conduction heat transfer in the crystal and deposition of the nutrient on the side 
walls. Details of the investigation appear in Mackowski (1996a) (Appendix B). The re- 
sults show that solid phase heat transfer, combined with the latent heat release from the 
depositing vapor, can significantly retard the overall growth rates in PVT systems. The 
predicted transport conditions under which heat transfer has a controlling role on crystal 
growth are consistent with a simplified, one-dimensional analysis of combined heat and 
mass transfer in PVT ampoules. The effect of heat transfer on the surface temperature 
distribution of the crystal is strongly dependent on the thermal boundary conditions at the 
side surfaces of the crystal substrate. Adiabatic conditions promote one-dimensional heat 
transfer, which leads to nearly isothermal surface temperature distributions and uniform 
growth rates. Cooled side walls, on the other hand, can lead to highly nonuniform surface 
temperature distributions, and correspondingly nonuniform growth rates. 
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Deposition of the crystal vapor onto the ampoule walls can have an equally signif- 
icant effect on the flux distribution. The redirection of vapor to the side walls has the 
primary effect of scavenging vapor from the crystal interface, which results in an overall 
decrease in mass transfer to the interface and a highly nonuniform flux distribution. For 
certain PVT conditions, our results also indicate that wall deposition can lead to buoyant 
recirculation even when the ampoule is oriented in a convectively ‘stable’ configuration. 
This phenomenon can result in interface flux distributions that are significantly different 
than those obtained under buoyancy-free, zero gravity environments for the same model 
conditions. 

5. Project Bibliography 

1. D. G. Walker, D. W. Mackowski, and R. W. Knight, “Effect of thermal creep on 
physical vapor transport crystal growth,” AIAA 32 nd Aerospace Sciences Meeting 
and Exhibit, 10-13 January 1994, Reno, NV, paper 94-0567. 

2. D. W. Mackowski and R. W. Knight, “Numerical Modelling of Physical Vapor Trans- 
port under Microgravity Conditions: Effect of Thermal Creep and Stress,” in Micro- 
gravity Studies of Organic and Polymeric Materials , D.O. Frazier and C.E. Moore, 
Eds., NASA Conference Publication 3250, pp. 17-20 (1994). 

3. D. W. Mackowski, R. V. Rao, and R. W. Knight, “Effect of solid phase heat trans- 
fer and wall deposition on crystal growth in physical vapor transport ampoules,” J. 
Crystal Growth , accepted for publication (1996). 

4. D. W. Mackowski, R. V. Rao, D. G. Walker, and R. W. Knight, “Numerical investi- 
gation of the effects of thermal creep in physical vapor transport,” J. Crystal Growth, 
submitted for publication (1996). 

6. Supported Personnel 

1. D. G. Walker (MS, 1994) 

2. V. R. Rao (Postdoctoral Research Associate) 


4 



Appendix A 


Numerical Investigation of the Effects of 
Thermal Creep in Physical Vapor Transport 

D.W. Mackowski, V.R. Rao, D.G. Walker and R.W. Knight 

Mechanical Engineering Department 
Auburn University, AL 36849 


Submitted to J. Crystal Growth, 2/96 


Abstract 

It has been recently recognized that the nonisothermal conditions present in physical 
vapor transport ampoules can give rise to a slip flow of gas over the side walls of the 
ampoule. This phenomenon, known as thermal creep, is usually insignificant relative to 
buoyancy-induced flows under similar nonisothermal conditions, and has therefore been 
neglected in previous PVT numerical models. However, thermal creep can, in principle, 
become a dominant convection mechanism in buoyancy-free environments such as those 
encountered in microgravity experiments. We present here a numerical investigation of 
the effects of thermal creep on the growth process in axisymmetric, binary component 
PVT systems. A nondimensional model, which includes buoyancy and Soret diffusion, is 
developed. We show that thermal creep can result in recirculating bulk flows within the 
ampoule. For relatively high values of the Schmidt number and large wall temperature 
gradients, these flows can result in significantly nonuniform distributions of mass flux at 
the crystal interface, and can also be comparable to or exceed the flow velocities generated 
by buoyancy under normal gravities. The effects of thermal creep on buoyant convection, 
and on the Soret transport of the vapor, are examined. 
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1. Introduction 


The formation of crystal structures is largely governed by the mass transport of the 
fluid material as it crystallizes at the nutrient/crystal interface [1]. This is certainly the 
case in physical vapor transport (PVT) crystal growth processes. In general, PVT is the 
process in which an amorphous solid is vaporized at one end of a closed, cylindrical am- 
poule, transported across the ampoule (usually through an inert carrier gas), and deposited 
onto a substrate that is maintained at a lower temperature than the source. The tempera- 
ture and/or composition gradients that are invariably present in PVT processes can result 
in vapor transport that is dominated by buoyancy-induced convection. However, under 
microgravity (/xg) environments the driving force behind buoyancy is largely eliminated 
and vapor transport occurs ostensibly in a diffusion-limited mode. Indeed, crystals grown 
in /xg PVT experiments have been shown to posess uniquely different and improved struc- 
tural properties that those grown on earth under similar operating conditions [2—4]. It is 
not known, however, to what extent the improvements result from modified vapor-phase 
transport conditions or the reduction in weight strain within the crystal. 

Because of the considerable cost in performing /xg PVT experiments, much effort has 
been given to numerical modelling of the process, Markham and Rosenberger [5] were 
the first to numerically examine the effects of microgravity on vapor transport rates in a 
cylindrical ampoule. Their results demonstrated that, under /xg conditions, the transport 
of vapor is not adequately described by one-dimensional diffusion of the nutrient through 
a stagnant carrier fluid. Rather, the rejection of the carrier at the source and crystal 
interfaces leads to a recirculation of the carrier within the ampoule, which in turn results 
in radial concentration gradients of the nutrient. The effect of carrier recirculation and 
buoyant recirculation on growth rates and growth uniformity have been further examined 
in comprehensive numerical investigations of PVT by Greenwell [6], Markham et al. [7], 
and Nadarajah et al. [8]. 

In general, a key conclusion of the previous modelling efforts is that the elimination 
of buoyancy-induced convection results in improved uniformity of growth at the crystal 
interface. Recently, however, Rosner [9] reported that the non-isothermal conditions en- 
countered in PVT processes can lead to additional, previously-unrecognized convection 
mechanisms. In particular, a temperature gradient tangential to a wall can result in a 
slip flow of gas over the wall. This phenomemon, which is known as thermal creep, is a 
macroscopic consequence of the molecular transfer of momentum to the wall. Molecules 
originating from the warmer regions in the gas will strike a wall with a higher momentum 
than those arriving from the cooler regions. The surface reacts to this uneven momentum 
transfer by imparting a force in the gas that is directed towards the warmer regions, which 
results in a slip flow over the surface. 
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The velocity of the creep flow is phenomenologically related to the temperature gra- 
dient by [9—11] 

c,v dT 

uc - V ' • ( 1 ) 


where v is the kinematic viscosity of the gas, z is the coordinate tangential to the wall, and 
C a is a dimensionless coefficient which depends on the momentum and thermal accommo- 
dation of the molecule-wall collisions. For diffuse reflection of the molecules the value of 
Cs is estimated as 2.2 [10]. Under normal gravity conditions thermal creep typically re- 
sults in a velocity that is significantly smaller than buoyancy-induced velocities produced 
under similar nonisothermal conditions, and has therefore been neglected in practically 
all numerical investigations of earthbound heat and mass transfer processes. In vanishing 
gravity, however, Rosner qualitatively indicated that the creep velocity can become com- 
parable to the typical Fickian diffusion velocity in PVT processes [9]. His conclusions were 
further supported by order-of-magnitude analyses by Napolitano et al. [12] and Viviani et 
al. [13]. 

To summarize these previous qualitative analyses, note first that the characteristic 
Stefan velocity of the vapor will be given by up = D abPz/T, where Dab is the binary 
diffusion coefficient, L is the ampoule length, and Pe is the mass transfer Peclet number 
based on the mass-averaged velocity. Assuming that the ampoule wall temperature varies 
linearly from the source temperature Ts to the crystal temperature Tc, the characteristic 
creep velocity u c at the crystal is obtained from Eq. (1) as uc « vC s {Ts — Tc)/LTc- 
Consequently, the ratio of the characteristic creep and Stefan velocities is 


Uc_ ~ uCa Ts - T C _ ScCa&£_ 
ud PcDab Tc Pe 


(2) 


in which Sc = v/ Dab is the Schmidt number and AT* is a characteristic nondimensional 
temperature difference. Considering the fact that AT* cannot realistically exceed values 
of around 0.5 (and will most likely be around 0.05-0.1), the creep and Stefan velocities 
will be comparable only for conditions of large Sc/ Pe. This condition will be attained 
for vanishing concentration of the nutrient in the carrier (Pe ■C 1) and, simultaneously, 
relatively large ratios of the molecular weights of the nutrient and carrier (large Sc). As 
noted in Refs. [9, 12, 13], there have been a number of /ig experiments that fall into this 
regime - most notably the copper phthalocyanine experiments of Debe et al. [3]. 

The objective of this work is to use detailed numerical modelling to quantitfy the 
effects of thermal creep on the PVT process. Since the conditions that are favorable to 
creep flows (i.e., large Sc and/or AT*) are also the conditions in which Soret (i.e., thermal) 
diffusion can become important, a second objective of the work is to examine the role of 
Soret diffusion (in combination with thermal creep) on vapor transport. Emphasis is given 
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to the effect of these transport mechanisms on the growth rates and growth unif or mi ty 
at the crystal interface, and the differences in the growth patterns that occur with and 
without buoyancy. 


2. Transport Model 

The geometry of the PVT ampoule is schematically illustrated in Fig. 1. The ampoule 
is cylindrical, with length L and radius R. Nutrient vapor, denoted as species A, originates 
from a planar source at z = 0 and deposits onto a planar crystal at z — L. The source and 
crystal interfaces are taken to be stationary, which is appropriate for the highly disparate 
vapor and solid phase densities of the nutrient. We also neglect the effects of heat transfer 
within the source and crystal, and take both surfaces to be at the uniform temperatures 
of Ts and Tc , respectively. A constant mass of carrier gas, denoted as species B, resides 
within the ampoule and is completely rejected at all surfaces. The system is axisymmetric, 
in steady-state, and all transport properties of the gas mixture are assumed constant. 

To facilitate the identification of regimes in which creep and Soret mechanisms con- 
tribute to vapor transport, we adopt a nondimensionai formulation of the governing equa- 
tions. The velocities are normalized by the characteristic diffusion velocity Dab/L, the 
temperature, pressure, and density of the gas are normalized by the values at the source, 
and the position variables are normalized with the ampoule length L. Denoting the nondi- 
mensionai variables with the superscript *, the equations for mass, momentum, energy, 
and species conservation become 


1 d(p*r*v*) d(p*u*) 

r* dr dz* 


= 0 


1 d(/>*r*u*u*) d(p*u*v*) 

r* dr* dz* 



+ Sc 


2 d , dv* 
r* dr* dr* 


2 d 

3 dr* 


V • v 


* 


(3) 


1 d{p*r*v*v*) d(p*u*u*) 

r* dr* dz* 


+ 


dz* 


( 


dv* du*' 
dz* dr* 


dP* 

-ATp^ + ScI 


1 d 


2u* 

(dv* 


dr \ dz * 


+ 



(4) 







+ p*g* 


1 gpWT) d(p*u*T*) 

r* dr* dz* 

1 d(p*r*v*w) d(p*u*w ) 

r * Q r * 4 " q z + 


Sc (1 d , &£_ d 2 T* \ 
Pr \r* dr* dr* dz* 2 ) 

-V • p*v* A . 


( 5 ) 

( 6 ) 
(7) 


In the above, u* and v* denote the z and r components of the nondimensionai mass 
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averaged velocity, and w is the mass fraction of the nutrient. The dimensionless parameters 
appearing in the momentum equations are defined 

PsL 2 . 9l 3 


Np = 


psUab 


0\b 


in which subscript S refers to the value at the source, and Pr = c p p,/k is the Prandtl 
number of the mixture. Note that the Schmidt number is defined with respect to the 
density at the source, i.e., Sc = p/ psP>AB, and the Rayleigh number for expansive/solutal 
buoyant convection is Roll = |(1 — pc/ps)d*\/Sc. The nondimensional diffusion velocity 
of the nutrient, v^, is comprised of Fickian and Soret components via 

vi = -V. - a T (9) 

in which a-r denotes the Soret (or thermal diffusion) factor [1,14]. The corresponding 
cross-product transport mechanism in the energy equation, i.e., Dufour diffusion (enthalpy 
diffusion via a concentration gradient), has been neglected - which will be appropriate 
for the Pr/Sc conditions of interest. Likewise, we have neglected the contribution of 
temperature gradients to fluid stress (‘thermal stress’) in the momentum equations. Based 
on the order-of-magnitude estimates in Refs. [12, 13], this mechanism will not significantly 
affect the flow field for the conditions examined here. 

We assume ideal behavior of the gas in the ampoule. We also neglect the effects 
of flow-induced pressure gradients on the gas density - which is entirely reasonable for 
the low Reynolds number conditions of typical PVT systems. The density will then be 
a function solely of temperature and composition. In nondimensional form, the ideal gas 
law will then become 

P ~ ~rjiP 1 ( 10 ) 

with a nondimensional molecular mass of the mixture, M * , given by 

M * M_ _ ws(l ~.Nm) + Nm 

Ms w(l ~N m ) + Nm' K > 

where Nm = Ma/Mb is the ratio of molecular masses of the nutrient and carrier. 

Infinite-rate deposition kinetics are assumed at the source and crystal interfaces, which 
will correspond to fixed mole fractions of the nutrient at the interfaces for a given total 
pressure. The boundary conditions at the source are then 


T*(r*,0) = 
w(r*,0) = 


1 + {Nm — 1 )x.s 


v*{r*, 0) = 0 

u*(r*,0) = — 


dw aTW dT* 
dz* T* dz* 
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in which xs is the specified mole fraction of the nutrient at the source. Analogously, the 
boundary conditions at the crystal interface are 


T*(r*,l) = ra 

(16) 

, , ,, Nmxc 

W(r ’ 1)= l + (N M - l)xc 

(17) 

v*(r * , 1) = 0 

(18) 

x 1 f dw arw dT*\ 

(19) 

U(r ’ 1)= 1 -w\dz‘ + T- dz- ),. =l ’ 


where = Tc/Xs. 

The temperature distribution on the wall plays a key role in the magnitude and result- 
ing effects of thermal creep. To provide a degree of freedom in choosing the temperature 
distribution - yet maintain a relatively simple form - we model the wall temperature in a 
piecewise-linear fashion. The distribution is described by two characteristic lengths; Zhot 
and z co id. For wall positions between z = 0 and Zhot the wall temperature is equal to the 
source temperature, between Zhot and z co u the wall temperature varies linearly from the 
source to the crystal temperature, and from z = z co id to L the wall temperature is equal 
to the crystal temperature. This distribution is somewhat representative of the three-zone 
arrangement employed in the NASA fig crystal growth furnace [15]. We also assume that 
deposition of the nutrient does not occur on the side wall. The boundary conditions on 
the side wall can then be stated in nondimensional form as 


T*(R*,z*)={ 


r i. o <**<*:., 

1 - (1 -Tg) 


* _ * 

z z hot * ^ ^ 

? z hot ^ ^ z cold 


zLu ~ z* 


cold * hot 

IIS, 4, m <**<i 


dw 


= 0 


r* = R‘ 
* 


dr* 


, .C.ScST * 
v'(R',z') = 0, 


r*=R‘ 


( 20 ) 


( 21 ) 

( 22 ) 

(23) 


in which R* = R/L. At the centerline of the ampoule we impose the usual symmetry 
conditions, which correspond to zero v* and zero radial gradients of u*, T* , and w. 

Solution of Eqs. (3-7) was accomplished using the SIMPLE algorithm as described by 
Patankar [16]. Discretization of the domain followed the staggered-grid procedure required 
by the SIMPLE method. A non-uniform grid was employed which concentrated grid points 
near the interfaces and walls of the ampoule. Typically, 23 grid points were used in the 
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radial direction, and up to 53 grid points were used in the axial direction. In all cases, 
the presented results for dimensionless velocities are within 0.5% of the values obtained 
when the number of grid points was doubled. Details of the solution method are given in 
Ref. [17]. 

3. Results and Discussion 

3.1. System parameters 

The relevant dimensionless parameters include those pertaining to the fixed boundary- 
conditions ( xs , xc, T£, zl ot and z* old ), those describing the thermophysical and transport 
properties of the mixture (5c, Pr, C s , &T and Nm), and those describing the ampoule 
geometry and environment ( R * and g*). The actual value of the pressure number Np is 
irrelevant (providing it is a positive number) since we assume pressure gradients have no 
effect on density. 

Fixed values of parameters used in this work are Pr = 1, C s — 2.2 (which represents 
the theoretical maximum), z* old — 1, and R* = 1/5. The Soret factor ap is set either 
to zero (no thermal diffusion) or to 35c/4 - which represents the limiting value attained 
for a dilute concentration of a high molecular weight nutrient in a low molecular weight 
carrier [14, 18]. 

It should be noted that, in principle, the ratio of nutrient mole fractions at the source 
and crystal interfaces could be related to the temperature ratio via a Clausius-Clapeyron 
equation for the specific crystal material and ampoule total pressure. We choose not 
to introduce such a relationship - since it would eliminate one parameter (e.g., xc) by 
introducing another (the enthalpy of vaporization of the crystal). We treat the mole 
fractions and the temperatures as independent - yet attempt to assign values to them that 
are in the realm of physical plausibility. 

3.2. Diffusive/thermal creep convective flows 

We begin by examining the effect of thermal creep on the flowfield and interface 
flux in the absence of gravity and Soret diffusion. As a starting point, we consider a 
configuration that is somewhat representative of the H 2 -I 2 system modelled in Refs. [5—8]. 
The molecular mass ratio is set to Nm = 254/2 = 127, and the source and crystal nutrient 
mole fractions are xs — 0.7 and xc = 0.3. The one-dimensional mass transfer Peclet 
number for this situation, which is defined [6] 



is 1.2. We take the dimensionless crystal temperature to be Tq = 350/400 = 0.875, which 
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represents a somewhat larger temperature difference across the ampoule than those used 
in the previous investigations. 

We focus first on the effects of an increasing creep velocity throughout the entire 
ampoule. To do this, we impose a linear wall temperature distribution from the source to 
the crystal (i.e., zl ot = 0 and z* old = 1) and vary Sc from 1 to 10 while holding all other 
parameters constant. This is, of course, somewhat inconsistent with physics - because 
an increasing Sc would generally be associated with an increasing Nm and a decreasing 
nutrient vapor concentration. Our rationale for this approach, however, is that the next 
level of approximation would require an accounting of the effects of mixture composition 
on transport properties. We wish to avoid over-specifying the system - and keep the 
model simple to the maximum extent possible. We also note that for the buoyancy-free 
and Soret-free systems, the dimensionless velocity results are not appreciably altered by 
the value of Nm- 

Figure 2 shows the velocity field that results from the combined actions of the Stefan 
flow and thermal creep at the side wall, with Sc = 1 on the left plot and 10 on the right. 
The scale of the vector arrows is the same in both plots. Note first that the case with 
Sc = 1 results in a velocity field that is practically indistinguishable from the zero-slip 
case. This is not surprising, because the creep/Stefan velocity ratio, defined in Eq. (2), 
is u c /ud « 0.2 for this case. Upon increasing Sc to 10, the creep velocity at the wall 
becomes comparable to the Stefan velocity. The combined effect results in a recirculating 
flow, with the cold fluid pumped towards the source along the walls, and returning to the 
crystal along the ampoule centerline. 

Even though thermal creep results in a roughly twofold increase in centerline velocity 
for Sc = 10, the total mass transfer rate to the crystal is only slightly increased (ca. 5%) 
by the creep-induced recirculation. The recirculating flow does, however, have a significant 
effect on the distribution of mass flux at the crystal interface. Shown in Fig. 3 are growth 
distributions as a function of radius for Sc ranging from 1 to 10, with all other model 
conditions the same as in the previous plot. We present a renormalized mass flux, which 
corresponds to the ratio of the dimensionless normal velocity at the interface divided 
by the 1-D Peclet number predicted from Eq. (24). For the Sc = 1 case the velocity 
distribution is essentially unchanged from that attained in the absence of creep. The fact 
that the distribution is not completely uniform in this case is due to the recirculation of 
the carrier in the ampoule, as discussed in Ref. [6]. An increasing magnitude of the creep 
velocity, brought on by large values of Sc, augments the radial segregation of the nutrient 
concentration above the interface. Centerline growth rates are seen to increase by roughly 
20 percent for the given range of Sc, which are balanced by a decrease in growth at the 
edges of the interface. 
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A linearly decreasing wall temperature distribution thoughout the entire ampoule 
is not representative of actual PVT conditions. Frequently, the wall is maintained at 
a temperature at or above the source temperature to a position just above the crystal 
interface - the point being to prevent the parasitic nucleation of the nutrient onto the 
wall [8]. With regard to thermal creep, this effect would increase the temperature gradient 
- and resulting creep velocity - at the wall adjacent to the interface. To investigate 
this effect, we changed the location of z^ ot to points progressively closer to the crystal 
while maintaining Sc — 5 and all other conditions the same as in the previous results. 
Plots of the velocity field for z\ ot — 0.75 and 0.875 are given in Fig. 4. A vortex is 
now created within the ampoule section corresponding to the wall temperature gradient, 
which increases in intensity as the length of the section decreases. Normalized growth 
distributions calculated for the changing wall temperature distributions are plotted in 
Fig. 5, with the same normalization used in Fig. 3. As can be seen, a steepening of the wall 
gradient leads to an increasingly nonuniform growth distribution - with the greatest effects 
occuring at the edges of the interface. The central factor here is that the creep-induced 
recirculation occurs directly above the interface - and thus has a direct influence on the 
transport of nutrient at the crystal-vapor boundary. This behavior could be significantly 
reduced or eliminated by minimizing the wall temperature gradient at the interface. 

S.S. Thermal creep and Soret diffusion 

The results of the previous section indicate that creep flows will have bearing on 
vapor transport only for relatively large Sc and small T£. Large Sc values are typically 
associated with a large nutrient /carrier molecular mass ratio and a vanishingly small mass 
concentration of the nutrient in the gas — which will yield the maximum values of the Soret 
factor ax- The ‘creep regime’ would thus be synonymous with a significant Soret transport 
of the vapor. This equivalence is not surprizing. Soret diffusion can be viewed as a 
molecular thermophoresis - especially in the limit of large nutrient molecules [14]. Thermal 
creep at the wall and Soret transport of the vapor originate from the same fundamental 
mechanism - i.e., temperature gradients - and conditions which are favorable to one will 
be favorable to the other. 

To model such conditions, we chose dimensionless parameters that were representative 
of the copper phthalocyanine PVT experiments conducted by Debe et al. [3]. Specifically, 
we used T£ = 360/600 = 0.6, xs — 10" 3 , xc = 10 -5 , Nm — 600/30 = 20, Sc = 10 and 
g* = 0. The one-dimensional Peelet number for this condition is Pe id ~ 0.02 Information 
on the wall temperature distribution is not obtainable from Ref. [3]; we therefore considered 
separate cases in which the location of z^ ot was moved from the source to the crystal. 

Similar to thermal creep, the effect of Soret diffusion on nutrient transport to the 
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crystal is strongly dependent on the thermal boundary conditions of the side wall. Unlike 
mass transfer, heat is transferred to (and, occasionally, from) the ampoule through the 
side walls. The temperature distribution in the gas - and its resulting effect on the con- 
centration distribution via Soret diffusion - can therefore be highly two dimensional. To 
illustrate this behavior, we show in Fig. 6 contour plots of the nutrient mass fraction w 
in the ampoule, calculated without thermal creep, for different locations of z^ ot . For the 
low-Pe conditions used here, the linear temperature distribution (zj£ ot = 0) is closely rep- 
resentative of an adiabatic wall. Heat transfer for this case is essentially one-dimensional, 
and acts to transport the nutrient uniformly to the crystal. Note also that the mass frac- 
tion field for the z£ of = 0 case is qualitatively similar to that attained under high Pe, 
convective/diffusive conditions - in that the nutrient tends to ‘pile up’ at the crystal in- 
terface. Indeed, the processes of convection (via the Stefan flow) and Soret transport of 
nutrient are analogous for the given conditions. Within the bulk gas nutrient transport 
is dominated by Soret diffusion. At the crystal, however, the contribution of Soret diffu- 
sion essentially vanishes (because of the relatively small mass fraction of nutrient at the 
crystal), and deposition becomes controlled by Fickian diffusion. 

As zl ot is moved away from the source and towards the crystal, a progressively larger 
fraction of the heat that is transferred to the crystal originates from the heated side walls of 
the ampoule. This creates a two-stage effect within the mass fraction field. For regions well 
within the heated wall zone the temperature is nearly uniform at Ts. Nutrient transport 
therefore becomes controlled within this zone by the relatively slow process of Fickian 
diffusion - which results in the high concentration gradients. At the transition point 
between the heated and gradient zones, radial heat transfer from the wall acts to drive the 
nutrient towards the centerline. At the same time, the large axial temperature gradients 
within the gradient zone accelerate the transport of nutrient, which results in the low 
concentration gradient within this second zone. 

The radial segregation of the nutrient that results from wall heat transfer can, as 
expected, have a significant effect on the flux distribution at the interface. This is demon- 
strated in Fig. 7, in which we plot the normalized interface velocities for different zj£ ot 
locations. Similar to the thermal creep effects discussed in Sec. 3.2., moving z\ oi closer to 
the interface results in progressively nonuniform, axial-peaked growth distributions. The 
extreme (and idealized) case of z* hot = 1 (i.e., a wall entirely at the source temperature) 
results in a complete absence of deposition for interface positions adjacent to the wall. 
This particular case actually resulted in physically-impossible results in the wall-interface 
corner (e.g., negative values of to) — which are a numerical artifact of the essentially in- 
finite temperature gradient at the corner. The effect of Soret diffusion on the total mass 
transfer rate, on the other hand, is greatest for z* hot = 0 (i.e., an adiabatic wall) - since this 
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condition results in an axial temperature gradient across the entire length of the ampoule. 

For the given model conditions, the thermal creep velocity is on the order of 500 times 
greater than the characteristic diffusion velocity. Consequently, including creep into the 
model results in convectively-dominated mass transport within the gradient section of the 
ampoule - as opposed to essentially convection-free conditions in the absence of creep. 
The recirculation generated by the creep has a controlling role on the distribution of flux 
at the interface, and the interfacial velocity profiles (shown in Fig. 8) attain a form that 
is qualitatively similar to those shown in Fig. 5. Soret diffusion, on the other hand, now 
affects primarily the overall mass transport to the interface. As was the case in the results 
of Sec. 3.2., inclusion of thermal creep has a relatively small effect on the total transfer of 
the nutrient. 

3. 4- Effect of creep on buoyant convection 

The large Sc and large AT conditions in which creep becomes significant would cor- 
respond to conditions in which gas density is controlled primarily by temperature. It is 
therefore relevant to compare creep to buoyant flows that axe driven by thermal expansion 
- as opposed to solutal-driven convection. To examine the interplay between buoyancy 
and creep, we performed calculations using the same conditions as in the previous section, 
except with T£ = 0.8 and z* hot = 0.75, and varied g* from 0 to — 5 x 10 7 . The negative 
value of g* indicates source-directed gravity, which represents the unstable configuration 
(for thermal expansion) of the crystal on top. As a point of reference, an ampoule of 
10 cm length, with Dab = 12.4 cm 2 /s (which corresponds to the CuPc value used in 
Ref. [3] at 5.4 torr total pressure) would have g* = —6.3 x 10 3 . Since Dab varies inversely 
proportional to total pressure, the chosen range of g* would represent operation of this 
experiment at total pressures up to about 0.6 atm. 

Before presenting results of our calculations, we should note that the axisymmetric 
formulation of our model is of questionable validity for buoyant convection in vertical 
cylinders [7]. However, in the absence of thermal creep, the heated side walls of the ampoule 
would be expected to stabilize an initial axisymmetric buoyant recirculation as g* decreases 
from zero. The gas, heated by the walls, would rise along the walls to the crystal and 
sink towards the source along the centerline. We have confirmed this behavior in limited 
calculations using a 2-D, cartesian representation of the ampoule (similar to that used in 
Ref. [8]), with equivalent interface and side wall boundary conditions. Once —g* exceeds 
a critical value the axisymmetric roll would become unstable, and would transform to a 
single, asymmetric roll. Based upon our 2-D cartesian model, the transformation to a single 
roll occurs for g* « —2 x 10 7 for the given model conditions. Since a cylindrical system 
would be expected to be more stabilizing than a 2-D cartesian system, the axisymmetric 
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assumption is probably valid throughout the given range of g* . Of course, validation of 
this assumption would require a 3-D model - which is beyond our current computational 
capabilities. 

A distinguishing feature between creep recirculation and expansive buoyant recircu- 
lation is that they occur in opposite directions. Creep drives the gas along the gradient 
section of the wall towards the source, whereas buoyancy is directed along the wall towards 
the crystal. The creep flow would thus be expected to inhibit the development of buoyant 
recirculation for this particular configuration. The interaction between creep and buoyant 
flows is graphically illustrated in Fig. 9, in which we present velocity vectors calculated for 
g* = — 5 X 10 5 (left), —5 x 10 6 (middle), and —2 x 10 7 (right). The scale of the velocity 
vector is the same in all plots. At the smallest value of g* the recirculation is completely 
dominated by the creep flow. At the next g* level creep still controls the convection within 
the gradient zone, yet buoyancy is seen to have created a counter-rotating vortex directly 
below the creep-produced recirculation. Finally, increasing g* to —2 x 10 7 results in an 
extremely complicated flow pattern, The upwards, buoyancy-controlled flow now extends 
to the crystal, and the creep-generated vortex has been pushed into a thin region adjacent 
to the side wall. 

The transition from creep to buoyancy controlled convection is reflected in the inter- 
facial velocity profiles, which are shown in Fig. 10 for the same conditions as in Fig. 9. 
Also shown are profiles calculated in the absence of thermal creep. The effect of combined 
creep and buoyant flows is seen to be qualitatively similar to that obtained by superim- 
posing the separate velocity profiles for pure creep and pure buoyant convection. This 
behavior is somewhat expected. Even though the creep flow has a major effect on nu- 
trient transport, the transport of energy and momentum are still controlled primarily by 
diffusion - for which a superposition of the independent creep and buoyant and flow fields 
becomes appropriate. This is further illustrated by examination of the total mass transfer 
rate to the crystal, denoted as j*, as a function of g*. We present this information in 
Fig. 11 for the zero-slip and thermal creep predictions, in which j* has been normalized 
by jl D — nR* 2 Pe lD . Thermal creep is seen to slightly delay the onset of buoyant con- 
vection control of mass transfer - which is indicated by the sudden increase in j* near 
g * « —2 x 10 7 - yet for the most part the trends with and without thermal creep are 
similar. 

Calculations were also performed, using the same conditions, of creep in conjunction 
with buoyancy stabilized (crystal on bottom) conditions. The results for this case show 
that the effects of creep are not appreciable altered by the presence of stabilized gravity 
relative to that calculated for zero gravity. 


A-12 



4. Summary 

The results of this investigation reveal that, for certain operating conditions, thermal 
creep at the ampoule side walls can significantly re-distribute the mass flux arriving on 
PVT crystal interfaces. In general, the investigation supports the conclusions drawn from 
the order-of-magnitude estimates of Refs. [12,13], in that the effects of creep become 
most evident for relatively large Sc, large AT and small Pe conditions. Indeed, flows 
resulting from normal-gravity buoyant convection can become, in this regime, relatively 
minor compared to the often-assumed ‘higher-order’ mechanism of thermal creep. 

The model we have presented, however, is highly idealized, and is not meant to be a 
precise prediction of growth processes in actual PVT experiments. The crystal interfaces 
in such experiments will not be isothermal because of latent heat release by the depositing 
vapor and finite conduction heat transfer resistance in the crystal [17, 19]. The interfaces 
will also not remain planar - especially in the presence of mechanisms (such as buoyancy, 
thermal creep and Soret diffusion) that tend to nonuniformly distribute the flux at the 
interface. Another uncertainty in the present analysis are the values of the slip coefficient 
on the side walls ( C a ) and the Soret coefficient {cut). Our model predictions used the theo- 
retical maximum values of these properties, and consequently the actual effects of thermal 
creep and Soret diffusion will likely be smaller than those predicted here. Nevertheless, 
our results strongly indicate that these frequently-disregarded mechanisms can have bear- 
ing on vapor transport - and should therefore be considered in detailed, comprehensive 
modelling efforts of actual PVT processes. 
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Figure Captions 


1. Ampoule Schematic 

2. Dimensionless velocity vectors for Sc = 1 (left) and Sc = 10 (right), for a linear tem- 
perature gradient across entire ampoule side wall. Other model parameters described 
in text. Origin of domain is lower left corner. 

3. Normalized interfacial velocities vs. dimensionless radial position r/R as a function of 
Sc, other conditions the same as in Fig. 2. 

4. Dimensionless velocity vectors for z^ ot = 0.75 (left) and = 0.875 (right), with 
Sc — 5 and all other parameters the same as in Fig. 2. 

5. Normalized interfacial velocities vs. dimensionless radial position as a function of Z hoV 
other conditions the same as in Fig. 2. 

6. Constant-value contours of nutrient mass fraction for z\ ot = 0 (left), z* hot — 0.5 
(middle), and = 0.875 (right). Zero slip boundary condtions, diffusive and Soret 
transport of nutrient. 

7. Normalized interfacial velocities vs. dimensionless radial position as a function of z^ ot 
for the conditions of Fig. 6. 

8. Normalized interfacial velocities vs. dimensionless radial position for the conditions of 
Fig. 7, with the addition of thermal creep convection. 

9. Dimensionless velocity vectors for g* = —5 x 10 5 (left), — 5 x 1.0® (middle), and —2 x 10 7 
(right). 

10. Normalized interfacial velocities vs. dimensionless radial position as a function of g* , 
for thermal creep and zero slip conditions. 

11. Normalized total mass transfer rate vs. g* for thermal creep and zero slip conditions. 
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Fig. 2, Mackowski et al. 
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Abstract 

A numerical investigation of binary component, convective/ diffusive physical vapor 
transport in cylindrical ampoules is presented. We specifically examine the effects of 1) con- 
duction heat transfer within the crystal, and 2) deposition of the vapor onto the ampoule 
side walls, on the growth rates and growth uniformity of the crystal. A two-dimensional 
numerical model of heat and mass transfer within the crystal and vapor regions of the 
ampoule is developed. Calculations are performed for H 2 -I 2 transport systems in unit and 
zero gravity conditions. Results indicate that the latent heat release from the deposit- 
ing vapor can significantly affect the interface temperature distribution of finite-thickness 
crystals. Depending on the thermal boundary conditions at the crystal side walls, heat 
transfer within the crystal can either increase or decrease the uniformity of mass flux at 
the interface. The predicted transport conditions under which heat transfer has a control- 
ling role on crystal growth are consistent with a simplified, one-dimensional analysis of 
combined heat and mass transfer in PVT ampoules. Calculation results also reveal that 
deposition of the vapor onto the side walls of the ampoule can lead to highly nonuniform 
distributions of mass flux at the crystal interface. Furthermore, side wall deposition can 
result in buoyant recirculation even when the ampoule is oriented in a convectively stable 
configuration. This phenomenon can result in interface flux distributions that are consid- 
erably different than those obtained under zero gravity environments for the same model 
conditions. 
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1. Introduction 

In the most general sense, crystal growth by physical vapor transport (PVT) involves 
the steps of 1) sublimation of a nutrient vapor from an amorphous source 2) transport 
of the vapor through an inert carrier gas to a substrate that is maintained at a lower 
temperature than the source, and 3) deposition of the vapor onto the substrate forming a 
crystal. As noted by Rosenberger [1], the quality of crystals that are grown in this process 
are governed largely by the heat and mass transfer mechanisms that occur within the PVT 
ampoule. Since the density of the gas can vary significantly within the ampoule, buoyancy- 
driven convection can often become a dominant mode of vapor transport. Numerical 
investigations by Markham et al. [2], and more recently Nadarajah et al. [3], have quantified 
the effects of buoyant convection on the flux distribution in PVT ampoules. As expected, 
buoyant recirculation of the vapor within the ampoule can lead to highly nonuniform 
distributions of mass flux at the crystal surface. 

To reduce this effect while maintaining relatively high transport rates, the potential 
of manufacturing crytals under microgravity (jig) conditions has received a fair amount 
of attention. Indeed, crystals grown in jig PVT experiments have been shown to posess 
uniquely different and improved structural properties than those grown on earth [4—6]. 
At this time, however, it is not certain whether the differences result from modified vapor- 
phase transport conditions, or are due to the elimination of weight stress in the crystals. 
The combined numerical investigations of Greenwell et al. [7] and Markham et al. [2] would 
indicate that vapor transport under convectively-stabilized, unit gravity conditions (i.e., 
when the high-density vapor is on the bottom) is not markedly different than zero-gravity 
conditions. 

The numerical investigations in Refs. [2, 3, 7], however, were limited to conditions of an 
isothermal crystal interface and zero deposition on the side walls of the ampoule. As noted 
in Ref. [3], an isothermal interface would not be realistic for poorly conducting crystals. 
Latent heat release of the depositing vapor can result in a large heat flux at the interface, 
and the interface temperature distribution would be dependent on the distribution of the 
heat flux as well as the removal of heat from the crystal. Since the saturation vapor pressure 
of the nutrient is often highly sensitive to temperature, a relatively small variation in 
temperature across the crystal surface could have an appreciable effect on the local growth 
rates. It has been observed that conduction of heat from the crystal to the side walls of 
the ampoule can lead to a preferential deposition of vapor at the wall, ultimately resulting 
in a concave interface shape [8]. 

The purpose of this paper is to examine the roles of heat conduction in the crystal and 
deposition onto the side walls on the growth process in convective/diffusive PVT systems. 
We begin by presenting a one dimensional analysis of combined heat and mass transfer in 
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PVT ampoules, which will provide a simple estimate of the conditions at which growth 
becomes controlled by heat conduction in the crystal. Next, we use a two-dimensional 
numerical model to predict the effects of crystal conduction and wall deposition on the 
mass flux distribution for unit and zero gravity environments. 

2. Model 


2.1. One dimensional heat/mass transfer 


An analysis of pure mass transfer in one-dimensional, steady state PVT systems 
has been presented in Ref. [7] . More recently, Louchev and Zaletin [9] developed a one- 
dimensional model for combined solid phase heat transfer and gas phase mass transfer 
in small Peclet number regimes (i.e., dominant diffusive transfer of the vapor). Here, we 
extend the one-dimensional model to include convective/diffusive heat and mass transport. 

The geometry of the PVT ampoule is schematically illustrated in Fig. 1. The ampoule 
is cylindrical, with length L and radius R. Nutrient vapor, denoted as species A, originates 
from the source at z = 0 and deposits onto the crystal at z = L. A constant mass of carrier 
gas resides within the ampoule, denoted as species B. The thickness of the crystal is S. 
Temperature and mole fraction are assumed to vary only in the 2 -direction, and all other 
properties are constant. We assume that the translation rate of the crystal interface is 
significantly smaller than the mass-averaged velocity of the vapor (which is appropriate 
for the highly disparate densities of the solid and gas) so that the system can be considered 
to be in steady state at any instant. The temperature of the source is Ts, and the ending 
wall is maintained at Tc,o- 

Unlike the model of Louchev and Zaletin [9], we restrict our analysis to infinite-rate 
deposition kinetics at the source and crystal interfaces. The nutrient vapor partial pressure 
at the interfaces will then be given by the equilibrium values at the interface temperatures. 
We also neglect radiation heat transfer within the ampoule. These restrictions will be 
appropriate for the carrier partial pressures (~ 100 torr), operating temperatures (~ 370 
K), and crystal materials (iodine) that are the subject of this investigation. 

The mass transfer Peclet number will be given by [7] 

u*L , / 1 - £a,« 


Pe = 


In 


Dab V 1 ~ x a,s ) ’ 

where u* is the mole-averaged velocity, Dab is the binary diffusion coefficient, xa is the 
mole fraction of the nutrient, and subscript i denotes the value at the crystal interface. 
The mole fraction at the interface can be related to the interface temperature Tc,i by 
evaluation of the Clausius Clayperon equation between Tc,,- and Ts; 

In (^iV=^ r (l~~ 

\ X A,S J RuTs \ Tc,i 


( 2 ) 
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in which h* g is the molar latent heat of sublimation for species A and R u is the universal 
gas constant. Combining the two previous equations results in an expression for Pe in 
terms of nondimensional thermal and thermodynamic properties of the system; 


Pe = In 


( 1 ~~ g A,s ex P IzEMU ~ ^«)] 

\ 1 — %A,s 



(3) 


in which H = h* g /R u Ts and 0,- = 1 — Tc,i/T$. 

The crystal interface temperature will depend on the dynamics of heat transfer through 
both the gas and the crystal. The gas and solid phase conductive fluxes at the crystal in- 
terface are related by 

9." ~ <lg - Cu*h sg . (4) 

The temperature distribution in the solid can be taken to be linear (1-D conduction), 
for which q " ■== k s (Tc,i — Tc,o)/S. The gas phase conductive flux can be obtained from a 
one-dimensional analysis of diffusive/convective heat transfer in the ampoule, which yields 

„ k g (Ts-T c ,i) PeLe 

q o L 'l-exp {-PeLe)' K) 


where Le = k g /Cc*D ab is the Lewis number, and C, k g and c* are the molar concentra- 
tion, thermal conductivity and molar specific heat, of the gas. With these substitutions, 
the nondimensional form of Eq. (4) becomes 

1 -+). ( 6 ) 

In the above, <j> = ( k s /S)/(k g /L ) is the ratio of gas and solid conduction heat transfer 
resistances, 9q = 1 — Tc,o/Ts , and 7 is the ratio of specific heats of the gas. Taken 
together, Eqs. (3) and (6) form a nonlinear relation for the Peclet number as a function 
of the nondimensional parameters <f>, 60, H, Le , 7, and xa,s- The mole fraction of A at 
the source can be eliminated using the Peclet number for Tc,i = Tc,o (i.e., zero crystal 
thickness), denoted as Pe 0, via 

= 1 ~ exp {-Pe 0) 

A ' S 1 - exp [— H6 0 /(l — 0 O ) — Peo] 

As <f> decreases, the mass transfer process will become limited by conduction heat 
transfer in the solid. In the limit of 4 > -4 0 (for which Oi -A- 0 and Pe 1), Eq. (6) predicts 
that Pe <x (j>, i.e., mass transfer will be proportional to k s and inversely proportional to 8. 
To identify the transition from mass to heat transfer control, we divide Eq. (6) by Pe, and 
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take the limit of <f> -4 0 while realizing that <f>/ Pe — > constant in this limit. After some 
manipulations, we arrive at 

+ Le(l-l/ 7 )tf)(i) 

1 ftcrit 

~ Pe 0 ' 

By plotting Pe/Pe o vs. 4>/4>crit we arrive at a set of similar curves for various thermo- 
dynamic and transport conditions of the system. These curves are presented in Fig. 2. 
The curves are based on thermodynamic property values for iodine (see Table 1), Lewis 
numbers of unity, Tc,o/Ts = 0.95, and Pe o of 0.01, 1, and 10. As expected, the curves 
asymptote to Pe/Pe o —¥ 1 and for <f>/<j> C rit > 1 and <§C 1, respectively, and the 

transition from mass to heat transfer control occurs for 4>/<f> C rit ~ 1. The behavior of the 
curves about 4>/4> C rit — 1 is controlled primarily by Pe o - in that H, do and Le have only 
a higher-order effect. For relatively small values of Pe o this transition is gradual, whereas 
for large Pe o sharp break occurs between the two regimes. 

One pertinent application of this model is a prediction of the crystal thickness at which 
growth goes from mass to heat transfer control. A characteristic value of this thickness 
is obtained by setting <f> = <f> C rit- Using the iodine property values given in Table 1 and 
Pe o = 1, this would correspond to S/L « 0.2. Actual crystal thickness - ampoule length 
ratios on this order would not be unrealistic. The inclusion of crystal heat conduction 
into a more detailed modelling of PVT systems for this particular type of crystal would 
therefore be warranted. 

2.2. Two dimensional heat/mass transfer model 

Similar to Refs. [2,3,7], we retain a primitive variable (i.e., pressure and velocities) 
formulation of the continuity and momentum equations. The computational domain is 
again illustrated in Fig. 1, The following assumptions are applied: 

1. The flow within the ampoule is axisymmetric. 

2. Gravity (if present) is directed towards the source. 

3. The carrier gas is completely rejected at the interfaces. 

4. Soret and Dufour effects are neglected. 

5. Transport properties (c p , k g , Dab , and the dynamic viscosity /i) are constant. 

Even though we seek steady-state solutions, we have adopted a transient numerical 

formulation of the governing differential equations. Our rationale is that the transient 
approach avoids many of the problems associated with the direct, iterative solution of 
the steady-state equations - such as insuring global mass conservation at each iteration 


( 8 ) 


A. - ( 1 ~ ex P [-ffflo/(l - flp )] 

Pe V #(exp (Peo)-l) 
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stage [3, 10]. With the given assumptions, the governing equations for continuity, momen- 
tum, energy, and species in the vapor phase region axe 
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( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 

(16) 
(17) 


In the above, p, P , and T are the density, total pressure, and temperature of the mixture, 
u and v are the z and r components of the mass-averaged velocity, and wb is the mass 
fraction of the carrier species. We pose the species equation in wb, rather than the common 
practice of wa = 1 — wb, for numerical reasons which will be described in the following 
section. We assume ideal behavior of the gas mixture, 


p _ pRuT 

M ’ 

in which the average molecular mass M is given by 


(18) 


M = 


1 ~ Wb WB 
M~A + Mb 


(19) 


The surface temperature of the source is taken to be uniform at Ts. The crystal 
surface temperature, on the other hand, is not specified a priori. Rather, we predict the 
crystal surface temperature by applying the conduction equation to the layer of crystal 
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material, again of thickness 5. Specifically, within the crystal the temperature (denoted as 
Tc ) is governed by 


dT c f (1 d dT c , d 2 T c \ 

iV> " dt ~ k ‘\rdr T dr + )' 


( 20 ) 


Assuming that the crystal/vapor interface is stationary and planar (i.e., not significantly 
deformed from a z = constant surface) the interface boundary condition is 


k 3 -— — k g — - puh 3 g, at 2 = L, (21) 

in which the latent enthalpy of species A, h 3g , is now on a unit mass basis. The temperature 
on the ending wall is taken to be uniform at the value Tc,o* 

As before, we set the partial pressure of the nutrient vapor at the crystal/vapor 
interfaces to the local saturation pressure. The boundary conditions at the source are then 
given by 


w B (r,0,t) = 1 - 


MA^A^satjTs) 


M b P 

u(r,0,t) = 0 

Dab dw B I 


1 - 


PAjSat^Ps) ( y _ Ma 

P V Mb 


-l 


u(r,0,t) = 


w b dz 


I z=0 


( 22 ) 

(23) 

(24) 


Completely analogous boundary conditions are obtained at the crystal interface. Note, 
however, that the saturation pressure P sa t at the interface will be a function of r - by 
virtue of the non-isothermal interface conditions. A specific relationship for the saturation 
pressure is presented in a following section. 

Temperature boundary conditions at the side walls will depend upon the specifics of 
the furnace system used to heat the ampoule. Here we assume that the wall temperature 
distribution follows that designed for the NASA /ig crystal growth furnace [11]. This system 
consists of three independent sections; a hot zone under which the wall temperature will 
be uniform at Ts; an adiabatic zone; and a cold zone with wall temperature of Tc,o- 
Consequently, the wall temperature boundary conditions are 


T(R,z,t) = T s , 0 < 2 < z hot 
= 0, Zhot Z <C. Zcold 

r—R 

T(R, 2 , t) = T c , o, z co i d <z<L + 8 , (25) 


dT 

dr 


in which Zhot and z co id are specified transition points between the furnace sections. 
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The boundary conditions for mass fraction and radial velocity at the side wall depend 
on the wall values of nutrient partial pressure and saturation pressure. For Pa/ Pa, sat < 1 
deposition does not occur and the boundary conditions are 


v(R,z,t) = 0 


8wb 

dr 


= 0. 

r=R 


(26) 

(27) 


If, on the other hand, the calculations predict that Pa/ Pa, sat > 1, we assume that depo- 
sition takes place. The boundary conditions will then become 


ws(R,z,t) 

v(R,z,t) 


MaP. A,3at[P(P; Z, t)] 

M B P 


Dab 

wb 


dwB 


dr 


r=R 


\ P A ,sat[T(R,Z,t )] 

(i Ma ^1 

L p 

V Mb) J 


(28) 

(29) 


The common no-slip boundary conditions (u(z, R, t) — 0) are applied at the wall, although 
it should be noted that wall temperature gradients can result in a thermal creep flow over 
the surface [10,12]. Calculations here will be limited to conditions in which creep has 
a negligible effect on vapor transport. Finally, symmetry conditions are applied at the 
ampoule centerline (r = 0), which correspond to zero v and zero radial gradients of u, T 
and t os- 

The system is initially taken to be at equilibrium at the crystal temperature Tc,o- 
Starting from t = 0, the temperature Ts is increased to the specified steady state value 
within a time period of t c , which was typically taken to be 0.1 s. The initial conditions 
and the transient heating process were chosen simply to fascilitate the calculation of the 
steady state condition in the ampoule. 

The final variable to be considered is the total pressure in the ampoule. Because 
velocities are specified at all surfaces, closure of the momentum and continuity equations 
alone do not require a unique specification of the pressure within the ampoule. That is, 
solution of these equations (for a fixed density field) will yield velocities u and v and a 
relative pressure distribution P'(r, z, t) in the ampoule. However, the total (or absolute) 
pressure distribution in the ampoule, at any point in time, is needed to calculate the 
density from Eq. (18). To determine the total pressure, we first adjust P' by a constant 
so that the spacial average of P' is zero. The total pressure is then obtained from P' plus 
the mean absolute pressure P m , 


P(r, 2 , t) = P m (t) + P'(r, z, t). 


(30) 
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The value of P m (t) can be determined from total mass conservation within the ampoule. 
Integration the continuity equation over the ampoule volume yields 


d_ 

dt 


n R pR pL 

prdrdz = {pu\ z=Q - pu\ z=L )rdr -R pv\ r=R dz. 
Jo Jo 


(31) 


The density can be eliminated using the ideal gas law, Eq. (18), in which P would be given 
by Eq. (30). The above equation will then yield a relation for the change in mean pressure 
P m as a function of u, v, P', T and wb- 


2.3. Numerical method 

We have employed the popular and well-established SIMPLE algorithm to calculate 
the flow and pressure fields. The discretization of the continuity and momentum equations 
follows the staggered-grid procedure required by the SIMPLE method [13], and the control 
volume approach was used to discretize the species equation and the gas and solid phase en- 
ergy equations. The interface conductivity method was used in the energy control volumes 
adjacent to the gas/crystal interface [13], with the latent heat contribution appearing as a 
source term. Time discretization was performed using a fully-implicit, forward-difference 
scheme. We use the discretized form of Eq. (31), combined with Eq. (18), to calculate the 
mean pressure P m at each time step, and then use Eqs. (30) and (18) to re-calculate the 
density within the ampoule. Note that this scheme provides a density field that identically 
satisfies global mass conservation for an arbitrary velocity field. This behavior is critical 
for the convergence of the SIMPLE algorithm. 

A subtle aspect in the numerical solution of the PVT system - which relates to our 
reason for using as the species variable - regards the proper application of the mass 
fraction boundary conditions at the solid/vapor interfaces. The mass fraction of the carrier 
is specified on the interfaces - yet the total flux of carrier (diffusive and convective) is also 
specified, i.e., zero. It is essential to maintain this latter condition during the iterative 
solution of the discretized species equation. Otherwise, the total mass of carrier will 
change in time - which is physically impossible. To meet this condition, it is necessary 
to discretize the transient and convective terms in the species equation in a divergent 
form (which corresponds to the form of Eq. (13)), as opposed to the conservative form 
recommended by Patankar [13]. A consequence of this approach is that the interfacial 
values of mass fraction no longer explicitly appear in the discretized species equation. 
Rather, they become implicitly coupled to the species equation solely through the velocity 
boundary conditions, i.e., Eq. (24). Aside from this modification, iterative solution of the 
equations follows the procedure outlined by Patankar [13]. The steady state condition 
is established when the change in the variables from one time step to the next becomes 
negligible. 
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One drawback of the transient formulation is that it is not possible to specify before- 
hand the steady-state mean pressure of the PVT system. Rather, the steady-state value 
of P m depends entirely on the specified initial and boundary conditions of the system. In 
situations in which the steady-state PVT behavior for a fixed P m is desired, it is thus nec- 
essary to guess an initial charging pressure P m , o that will yield the sought final pressure. 
As expected, this process involves iteration between P m ,o and the calculated P m . 

To verify our model we performed calculations to reproduce previous benchmark nu- 
merical investigations of two dimensional forced and free convective heat transfer in closed 
cavities [14,15]. Further validation of our model was obtained by comparing results for 
PVT velocity fields and growth rates to those calculated in Refs. [2, 3,7]. The calculated 
velocities and growth rates were found to be in acceptable agreement (ca. 5%) with the 
previous works. 

3. Results and Discussion 

3.1. System parameters 

To be consistent with the previous modelling efforts in Refs. [2, 3, 7], we focus our 
attention on the H 2 -I 2 PVT system. Values of the parameters in the PVT model that axe 
fixed thoughout this investigation are given in Table 1. The equilibrium partial pressure 
of I 2 vapor was given by the relation [16] 

Pa, sat = 101.32 x 10 3 exp ^30.795 - — - 2.013 In Pa. (32) 

From the Clausius-Clayperon equation, the latent enthalpy of the vapor at a characteristic 
interface temperature of 360 K is calculated from the above as h sg = 2.64 x 10 s J/kg, 
which agrees well with the value of 2.46 x 10 5 J/kg used by Abernathy et al. [8]. The initial 
charging pressure of the system (i.e., the pressure when the system is at equilibrium at 
Tc, 0 ) was set to 10 kPa. This resulted in a final steady-state pressure within 10% of 13.3 
kPa (100 torr) for all conditions examined. 

Using these parameters and the conditions of zero gravity and isothermal interfaces, 
a grid dependence study revealed that the differences between the calculated steady-state 
growth rates for uniform 12x12 and 22x22 grids (in the gas phase region) were normally 
within 3 percent, and the difference between 22x22 and 42x42 grids was within 1 percent. 
All the results presented here were calculated with the 22x22 grid. Depending on the 
thickness of the crystal, between 5 and 10 rows of grid points were used in the discretization 
of the solid-phase energy equation. 
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3.2. Effect of conduction in crystal 

Results of our calculations for combined solid and gas phase transport appear in Figs. 3 
and 4, in which we present distributions of steady-state mass flux at the crystal surface for 
a solid thermal conductivity of k so i = 0.5 W/m-K (solid iodine), with the crystal thickness 
5 randing from 0 to 5 mm. Results are given for zero gravity and unity gravity pointing 
towards the source (convectively stabilized). Figure 3 has a fixed side wall temperature 
of Tc, o in the crystal region, whereas the crystal side walls for Fig. 4 are adiabatic. For 
all cases the temperature of the side walls in the gas-phase region (0 < z < L) was 7s, 
and the length L of the vapor region was fixed at 2 cm. A constant value of L would be 
expected since growth of the crystal layer would be balanced by regression of the source 
layer. For the zero crystal thickness, zero gravity case one sees the behavior identified by 
Greenwell et al. [7], in that flux is highest at the centerline, whereas stabilizing gravity 
results in a more uniform distribution for an isothermal interface as reported by Markham 
et al. [2]. 

A quick inspection of the curves in Figs. 3 and 4 reveals that the thermal boundary 
conditions on the side walls of the crystal can have a profound effect on the flux distribution. 
The constant temperature case (Fig. 3) is discussed first. As 8 increases, the flux for both 
gravity conditions decreases at the centerline and dramatically increases at the wall. The 
behavior at the centerline follows from the one dimensional analysis in Sec. 2.1. Regions 
of the crystal near the centerline are cooled primarily by conduction to the end wall. 
Increasing the crystal thickness acts to increase the conduction resistance, and therefore 
increases the crystal surface temperature. The behavior near the wall, on the other hand, 
is somewhat counterintuitive. Regions of the crystal near the wall are thermally coupled 
to the side wall, and therefore maintain temperatures near 7c, o regardless of the crystal 
thickness. Because of this, one might expect the flux at the wall to remain near the 
value predicted for zero crystal thickness - yet the flux actually increases with thickness. 
Transport mechanisms in the vapor phase are responsible for this behavior. The retardation 
of flux at the centerline acts to increase the nutrient concentration over the entire crystal, 
which in turn results in a higher concentration gradient - and mass flux - at the surfaces 
near the wall. 

The flux distributions for zero and unit gravity are nearly identical at the maximum 
crystal thickness of 5 mm. Growth rates are now controlled primarily by solid phase 
conduction, and differences in the vapor-phase transport conditions (i.e., buoyancy or the 
lack thereof) will therefore have minimal effect of the flux distribution. When the crystal 
side wall is adiabatic (Fig. 4) the effect of wall cooling is removed, and the flux profiles 
attain a nearly uniform distribution. As before, increasing crystal thickness lessens the 
differences in flux distribution between the zero and unit gravity conditions. 
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The calculated flux profiles shown in Figs. 3 and 4 are essentially unaffected by the 
thermal boundary conditions at the side walls of the gas-phase region. For the given 
transport conditions the heat flux to the crystal interface is dominated by latent heat 
release from the depositing vapor. Gas-phase heat conduction from the heated side walls 
to the interface will therefore have a minimal effect on the interface temperature. 

To illustrate the degree to which the interface can become nonisothermal, we plot in 
Fig. 5 interface temperatures as a function of radius for the same conditions as in Fig. 3. 
The variation in interface temperature for the 8 — 5 mm thickness is around 8 K across 
the 1 cm radius of the surface. Extremely high temperature gradients are predicted near 
the wall, which is a consequence of the large heat flux that results from latent heat release. 
Indeed, it seems unlikely that such conditions would be attained in practice, because the 
finite thermal conductivity of the wall material would result in wall temperatures somewhat 
elevated above the nominal values used here. Nevertheless, it is clear that the high heat 
fluxes resulting from vapor deposition can lead to significant temperature gradients within 
the crystal - especially when the side walls of the crystal are cooled. 

The effect of crystal conduction heat transfer on the overall mass transfer rates is 
presented in Fig. 6. Here we adopt the normalization used in the 1-D analysis, and plot 
the ratio Pe/Pe \-d,q vs. 4>/4> C rit, where Pe is the Peclet number calculated from the 2-D 
results and Pe\-D,o is that corresponding to the 1-D model with zero crystal thickness. 
The results are based on the isothermal crystal side wall, zero gravity configuration. Two 
sets of Peclet number ratios are presented, which correspond to Pe based on the mass aver- 
aged velocity and the centerline velocity. Also included are the 1-D analytical predicitions 
for the same thermal and transport parameters of the system. 

The numerical results for Pe based on the centerline velocity are seen to correspond 
reasonably well with the 1-D preditions over the entire range of crystal thicknesses. This 
is to be expected, because the interface temperature at the centerline is, again, controlled 
primarily by the axial conduction of heat across the crystal. On the other hand, the 
increased flux at the crystal edges - which results from the inherently two-dimensional 
nature of energy transfer in the ampoule - leads to mass-averaged Pe values that are only 
slightly affected by increasing crystal thickness. 

It is apparent from the preceeding results that heat transfer within a poorly conducting 
crystal can significantly affect the distribution of vapor deposition. The thermal boundary 
conditions on the side wall containing the crystal play a critical role in governing this 
distribution. Ideally, a uniform distribution would be sought, and this would be best 
realized by maintaining adiabatic side walls. This condition, however, is difficult to attain 
in practice [3]. Cooled side walls - which are more physically realistic - will eventually 
lead to a concave interface shape due to the preferential deposition of vapor at the crystal 
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edges. This has been observed experimentally by Abernathy et al. in the growth of iodine 
crystals in high Pe conditions [8], 

3.3. Effect of wall deposition 

To investigate the effects of wall deposition, we simply extended the cooled section of 
the side wall into the gas phase region. Results of the calculations appear in Fig. 7 and 8, 
in which mass flux is plotted vs. position for vapor-phase cold zone lengths ranging from 
A z co id — 2 to 10 mm. Figure 7 shows the flux at the interface (i.e., the intended target 
of deposition), whereas Fig. 8 shows the side wall flux. To focus solely on the effects of 
wall deposition, the end wall interface in both plots has been maintained at the uniform 
temperature of Tc,o- As before, results are given for zero and unit, stabilizing gravity. 

The first point to draw from Fig. 7 is that wall deposition results in a significant 
decrease in overall mass transfer to the crystal interface. The side wall deposition surfaces, 
being closer to the source, act to scavenge nutrient vapor from the crystal. To quantify this 
effect, the 2 mm section of wall deposition surface results in an approximately 30 and 50% 
decrease in the total mass flux reaching the interface for the zero and unit gravity cases, 
respectively. Equally significant is the effect of side wall deposition on the distribution of 
flux at the interface. A 90 percent variation in flux across the radius of the interface occurs 
for A z co id = 2 mm, as opposed to around 10 percent variation in the absence of side wall 
deposition. Not surprizingly, the location of maximum flux in the entire ampoule occurs 
on the cooled side wall positions that are nearest to the source (Fig. 8) - simply because 
they are exposed to a relatively higher nutrient concentration. 

Perhaps the most interesting aspect of side wall deposition axe its effects on buoyant 
convection in stabilized configurations. Our calculations indicate that side wall deposition, 
in conjunction with source-directed gravity, can act to create a recirculation flow within 
the gas. The nature of the recirculation is shown in Fig. 9, in which we present vector plots 
of the mass averaged velocity for unit, source-directed gravity. To explain the behavior of 
the flow, we note, first, that the vapor density at the wall deposition surface will be lower 
than that at the immediately adjacent centerline - assuming, of course, that Ma > Mb. 
Gravity will therefore impart a relatively higher body force in the gas at the centerline. 
Concurrently, a radial convective flow will occur as the vapor is swept to the side wall 
deposition surfaces - and the magnitude of this flow will be greatest at the boundary 
between the cold and hot zones of the ampoule. Once the length of wall exceeds a critical 
value (around 5-6 mm for the ampoule conditions used here) the combined effects of radial 
density gradients and radial velocities result in a recirculating flow. The recirculation is 
counter to that attained in destabilizing gravity (i.e., g pointed towards the crystal), in 
that the flow is directed towards the source at the centerline and towards the crystal at 
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the wall. Referring again to Figs. 7 and 8, the recirculation acts to significantly lower the 
interface flux at the centerline, and increase the flux at the hot/cold junction of the wall, 
relative to the zero gravity condition. 

We emphasise again that the recirculation shown in Fig. 9 results from a complicated 
interaction between the axial and radial-directed convective flow of the vapor and buoyancy 
forces. For example, regardless of the wall thermal conditions, removal of the Stefan flow 
(by setting normal velocities at all surfaces to zero) while maintaining all other conditions 
the same resulted in convectively stable systems for source-directed gravity. This is in 
contrast to unstable (crystal-directed) gravity, in which solutal density differences, by 
themselves, can result in convective recirculation within the ampoule [17] . 

5. Conclusions 

Our calculations certainly indicate that an accurate modelling of PVT ampoules 
should include heat transfer in the crystal and (if present) deposition on the side walls. 
Both of these factors have been shown to significantly affect the rate and uniformity of 
crystal growth. We also submit that the results could have bearing on the effects of buoy- 
ant convection on the morphological stability of crystals grown in PVT ampoules. Again, 
significantly different distributions in flux between zero and unit-stablized gravity were 
observed in the presence of side wall deposition. A primary factor behind these differences 
is the fact that the growth interface, as a whole, becomes non-planar. Because of this, it 
is conceviable that a non-planar crystal interface would evolve differently under stablized 
gravity and zero gravity conditions - regardless of the presence of side wall deposition. An 
accurate prediction of such effects would require a numerical scheme that can conform to 
the non-planar boundary of the crystal. We are currently engaged in the modification of 
our model to account for such effects. 
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Table 1. Model parameters 


M a 

254 g/mole 

Mb 

2 g/mole 

Dab 

1.25 cm 2 /s (@ 100 Torr) 

V 

1.5 xlO -4 g/cm-s 

Pr (= c P n/k g ) 

1 

7 

1.4 

c p 

3.5i? u 

T s 

378.1 K 

Tc, o 

358.1 K 

R 

1 cm 

L 

2 cm 

fag 

0,5 W/m-K 

hgg 

2.64 x10 s J/kg 
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Figure Captions 


1. Ampoule schematic. 

2. Ratio of 1-D Peclet numbers with and without solid phase heat transfer vs. dimen- 
sionless solid phase resistance parameter, 

3. Local mass flux at crystal surface vs. radial position, isothermal crystal side wall. 

4. Same as Fig. 3, except with adiabatic crystal side wall. 

5. Interfacial surface temperature distribution for conditions of Fig. 3. 

6. Ratio of Peclet numbers with and without solid phase heat transfer vs. dimensionless 
solid phase resistance parameter. 

7. Local mass flux at crystal interface vs. radial position. 

8. Local mass flux at side wall vs. side wall position. Conditions same as in Fig. 7. 

9. Velocity vectors for case of Figs. 7 and 8, unit stabilizing gravity. 
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Fig. 1, Mackowski et al. 
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Fig. 5, Mackowski et al. 










